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ABSTRACT 

Numerous cosmological hydrodynamic studies have addressed the formation of galaxies. Here we 
choose to study the first stages of galaxy formation, including non-equilibrium atomic primordial 
gas cooling, gravity and hydrodynamics. Using initial conditions appropriate for the concordance 
cosmological model of structure formation, we perform two adaptive mesh refinement simulations of 
^IO^A/q galaxies at high redshift. The calculations resolve the Jeans length at all times with more 
than 16 cells and capture over 14 orders of magnitude in length scales. In both cases, the dense, 
10^ solar mass, one parsec central regions are found to contract rapidly and have turbulent Mach 
numbers up to 4. Despite the ever decreasing Jeans length of the isothermal gas, we only find one site 
of fragmentation during the collapse. However, rotational secular bar instabilities transport angular 
momentum outwards in the central parsec as the gas continues to collapse and lead to multiple nested 
unstable fragments with decreasing masses down to sub- Jupiter mass scales. Although these numerical 
experiments neglect star formation and feedback, they clearly highlight the physics of turbulence in 
gravitationally collapsing gas. The angular momentum segregation seen in our calculations plays an 
important role in theories that form supermassive black holes from gaseous collapse. 
Subject headings: cosmology: theory — galaxies: formation — black holes: formation — secular 
instability 



L MOTIVATION & PREVIOUS WORK 

Since the first investigations of galaxy interactions 
(|Holmberglfl94l using light bulbs, the use of numeri- 
cal simulations in galaxy formation has developed dra- 
matically. Not only gravity but also hydrodynamics 
and cooling are standard ingredients in the sophisticated 
computer models studying galaxy formation and inter- 
actions. In hierarchical structure formation, dark mat- 
ter (DM) halos merge to form lar ger halos while the gas 
infalls into these po te ntial wells ([Peebles fc Dickd 119681 : 
IWhite fc Reeslll978D . IWhite fc ReesI provided~the basis 
for modern galaxy formation, in which small galaxies 
form early and continuously merge into larger systems. 

As more high redshift galaxies were observed in the fol- 
lowing 10 years, IWhite fc Fren^ (|1991h refined the theory 
to address the observed characteristics in these galaxies. 
In their model, the halo accumulates mass until the gas 
cools faster than a Hubble time, tn, which usually occurs 
when atomic hydrogen line, specifically Lya, cooling is 
efficient. This happens when the halo has Tvir > 10** 
K, where the cooling function sharply rises by several 
orders of magnitude because the number of free elec- 
trons able to excite h ydro gen greatly increases at this 
temperature (Spitzen il978t) . One can define a cooling 
radius, TcooI, in which the interior material is able to 
cool within a Hubble time. Once the halo reaches this 
first milestone, TcooI increases through additional accre- 
tion and cooling. A rapid baryonic collapse ensues when 
icooi ^ tdyn (Rccs & Ostrikcr 1977). The material ac- 
celerates towards the center, and its density quickly in- 
creases. In the model discussed in White & Frenk, this 
collapse will halt when one of the following circumstances 



occurs. First, angular momentum can prevent the gas 
from collapsing further, and the system becomes rota- 
tionally supported. Afterwards, this disc fragments and 
star formation follows. Alternatively, star formation does 
not necessarily develop in a disc component, but the en- 
ergy released by stars during their main sequence and 
associated supernovae (SNe) terminates the collapse. 

These concepts have b een applie d also to the earliest 
galax ies in the universe (iMo et all [1998: Oh & Haima^ 
2002t IBegelman et al.l 120061: Eodato fc NataraianI 
200d). Many studies (e.g. iQstriker fc Gnedin 19961 : 
Hai man et aLlI]^ ICenll200a ISomerville fc Livid f20?)l 
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(Wise fc AbellboOSD demonstrated that OB-stars within 
protogalaxies at z > 6 can produce the majority of 
photons required for reionization. These protogalaxies 
contain an ample gas reservoir for widespread star 
formation, and the accompanying radiation propagates 
into and ionizes the surrounding neutral intergalactic 
medium. Several high redshift starburst galaxies have 
been observed that host ubiq uitous star formation 
at z > 6 (iStanwav et~an l2003t fMobasher et"al] l2005l : 
iBouwens et al.ll2006[ ). Additionally, supermassive black 
holes (SMBH) more massive than lO^Mf?) are present at 
these redshifts (e.g. iBecker et all 120011: iFan et al.|[200l 
l2f)f)6h . Finally, a reionization signature in the polariza- 
tion o f the cosmic mic rowave background (CMB) at z 
~ 10 ([Page et al.ll2007t ) further supports and constrains 
stellar and SMBH activity at high redshifts. 

The distinction between SMBH formation and a star- 
burst galaxy should depend on the initial ingredients (i.e. 
seed BHs, metallicity, merger histories) of the host halo, 
but the evolution of various initial states is debatable. 
It is essential to study the hydrodynamics of high red- 
shift halo collapses because the initial luminous object(s) 
that emerges will dynamically and thermally alter its sur- 
roundings. For example, as the object emits ultraviolet 
radiation, the nearby gas heats and thus the characteris- 
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TABLE 1 

Simulation Parameters 



Name 


1 


^part 


grid 








Ax 




[Mpc] 












[R©] 


A 


1.0 


2.22 X 10^ 


44712 


1.23 


X 10* (498^) 


41 


9.3 X 10-3 


B 


1.5 


1.26 X lO'^ 


22179 


7.40 


X 10^ (420^) 


41 


1.4 X 10-2 



Note. — Col. (1): Simulation name. Col. (2): Comoving size of the simu- 
lation. Col. (3): Number of dark matter particles. Col. (4): Number of AMR 
grids. Col. (5): Maximum number of unique grid cells. Col. (6): Maximum level 
of refinement reached in the simulation. Col. (7): Resolution at the maximum 
refinement level. 



TABLE 2 
Halos of interest'' 



1 


z 


Mtot 


a 


Pc 




Tc 


Mbe 


[Mpc] 




[Mq] 




[cm~ 


'] 


[K] 


[Me] 


1.0 


15.87 


3.47 X lO'^ 


2.45 


5.84 X 


1021 


8190 


4.74 X 10^ 


1.5 


16.80 


3.50 X lO'^ 


2.59 


7.58 X 


1021 


8270 


1.01 X 10^ 



Note. — Col. (1): Box size of the simulation. Col. (2): Final redshift 
of simulation. Col. (3): Total halo mass. Col. (4): a of the total mass 
compared to matter fluctuations. Col. (5): Central halo density. Col. 
(6): Central gas temperature. Col. (7): Gravitationally unstable central 
mass. 

^ The subscript "c" denotes central quantities. 



tic Jeans mass increases, which may in hibit the accretion 
of new gas for future star formation (|Efstathioul 119921 : 
iThoul fc W einberg 1996). 

The following work will attempt to clarify early galaxy 
formation by focusing on protogalactic (Tvir > lO** 
K) halos and follow ing their initial gaseous collapse. 
IWise fc Abeil (|2007aL hereafter Paper I) studied the viri- 
alization of protogalactic halos and the virial generation 
of supersonic turbulence. In this paper, we address the 
gas dynamics of the continued, turbulent collapse of a 
halo and study the evolution and characteristics of the 
central regions. In later studies, we will introduce the 
effects from primordial star formation and feedback and 
H2 cooling. The progressive introduction of new pro- 
cesses is essential to understand the relevance of each 
mechanism. We argue that our results may be rele- 
vant for scenarios that envisage SMBH formation from 
g aseous colla pses. 

iLoeb fc Rasiol (|1994D and lBromm fc Loebl (|2003l) con- 
ducted smoothed particle hydrodynamics (SPH) simu- 
lations that focused on the collapse of idealized, iso- 
lated protogalactic halos. The former group concluded 
that a central IO^Mq SMBH must exist to stabi- 
lize th e thin gaseous d isc that forms in their calcula- 
tions. iBromm fc Loebl considered cases with and with- 
out H2 chemistry and a background UV radiation field. 
They observed the formation of a dense object with a 
mass M ^ IO^Mq, or >, 10% of the baryonic matter, 
in simulations with no or strongly suppressed H2 forma- 
tion. These calculations without metal cooling and stel- 
lar feedback are useful to explore the hy drodynamics of 
the collapse under simplified conditions. ISpaans fc Silkl 
([2OO6) analytically studied the collapse of 10'* K halos 
with an atomic equation of state. They find that ~0.1% 
of the baryonic mass results in a pre-galacti c BH with a 
mass ~1O^M0. iLodato fc NataraianI ()2006D also found 
that ^5% of the gas mass in M = IO'^Mq halos at 



z ~ 10 becomes un stable in a gaseous disc and forms 
a SMBH. Recently, IClark et all (|2007D studied the ef- 
fects of metal and dust cooling on the fragmentation of 
a collapsing protogalactic core with varying metallicities 
{Z — 0, 10"^, 10-^^0) and found the gas fragmenting 
ten times as much in the 1O"^^0 case than the primor- 
dial case. In addition, the fragments in the primordial 
case are biased toward larger masses. 

A runaway gaseous collapse requires angular momen- 
tum transport so material can inflow to small scales 
and form a central object. The stability of rotating 
gaseous clouds have been subject of much interest over 
the last four decades and was th oroughly detailed by the 
work of IChandrasekhad (|1969! . hereafter EFE) . In the 
1960's and 19 70's, studies u t ihzing virial tensor tech- 
niques fEFE: iLebwitg Il967l: lOstriker fc Tassoull fl96l 
'Ostriker fc Boidenheimer' '1973'), variational techniques 
(Lynden-Bcll fc Ostriker 1967;^Bardeen et al. 1977), and 
N-body simulations (jOstriker fc Peebles! Il973ij ) all fo- 
cused on criteria in which a stellar or gaseous system 
becomes secularly or dynamically unstable. The first 
instability encountered is an m = 2 bar-like instabil- 
ity that is conducive for angular mome ntum transport 
in ord er to form a dense, central object. Beg elman et al.l 
()2006l ) investigated the conditions where a gaseous disc in 
a pre-galactic halo would become rotationally ti nstable to 
bar formation (see IChristodoulou et al.l [l995a"g). They 
adap t the "bars within bars" scenario dShl osman et al.l 
1989, 1990), which was originally formulated to drive 
SMBH accretion from a gaseous bar that forms within 
a stellar galactic bar, to the scenario of pre-galactic BH 
formation. Here a cascade of bars form and transport 
angular momentum outwards, and the system can col- 
lapse to small scales to form a quasistar with runaway 
neutrino cooling, resulting in a central SMBH. The sim- 
ulations detailed below show how many central bar-like 
instabilities form. 
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In §2, we describe our siraulations and their cosmo- 
logical context. In the following section, we present our 
analysis of the halo collapse simulations and investigate 
the structural and hydrodynamical evolution, the initial 
halo collapse, rotational instabilities, and the importance 
of turbulence. In §4, we discuss the relevance of angu- 
lar momentum transport and rotational instabilities in 
early galaxy and SMBH formation. There we also exam- 
ine the applicability and limitations of our results and 
desired improvements for future simulations. Finally we 
conclude in the last section. 

2. SIMULATION TECHNIQUES 

To investigate protogalactic halo collapses in the early 
universe, we utilize an Eulerian structured, adaptive 
mesh refinen ient (AMR), cosnrological hydrodynamical 
code Enzo^ (jBrvan fc Normanlll997lll999i : IO'Shea et all 
|2004[ ). Enzo solves the hydrodynamical equations us- 
ing a second order accurate piecewise paraboli c method 
(|Woodward fc Colellal [l98l iBrvan et allll99l . while a 
Riemann solver ensures accurate shock capturing with 
minimal viscosity. Additionally Enzo uses a particle- 
mesh N-body method to calculat e the dynamics o f the 
collisionless dark matter particles (jCouchmanll 1 99 Ih . Re- 
gions of the simulation grid are refined by a factor of two 
when one or more of the following conditions are met: 

(1) Baryon density is greater than 3 times QbPoN^'^^~^'^\ 

(2) DM density is greater than 3 times ricDMPo^^'^^+'^', 
and (3) the local Jeans length is less than 16 cell widths. 
Here iV = 2 is the refinement factor; I is the AMR refine- 
ment level; — —0.3 causes more frequent refinement 
with increasing AMR levels, i.e. super-Lagrangian be- 
havior; po = SHq/SttG is the critical density; and the 
Jeans length, Lj ~ /^-KpG^LniH, where Hq, fc, T, 
p, p, and rriH are the Hubble constant, Boltzmann con- 
stant, temperature, gas density, mean molecular weight 
in units of the proton mass, and hydrogen mass, respec- 
tively. The Jeans length refinement insures that we meet 
the Truelove criterion, which requires the Jeans length to 
be re solved by at least 4 cells on each axis () Truelove et al.l 
Il997( ). Runs with a refinement criterion of 4, 8, and 16 
Jeans lengths have indistinguishable mass weighted ra- 
dial profiles. 

We conduct the simulations within the concordance 
ACDM model with WMAP 1 year parameters oi h = 
0.72, = 0.73, Qm = 0.27, = 0.024/i-2, and a 
primordial scale invariant (n = 1) power spectrum with 
0-8 = 0.9 (jSpergel et all l2003f ). h is the Hubble pa- 
rameter in units of 100 km s~^ Mpc~^. Ha, ^m, and 
fib are the fractions of critical energy density of vac- 
uum energy, total matter, and baryons, respectively, erg 
is the rms of the density fluctuations inside a sphere 
of radius 8h~^ Mpc. Using the WMAPl parameters 
versus the significan tly different WMAP third year pa- 
rameters (WMAP3; ISpergel et all |2007| ) have no effect 
on the evolution of individual halos that are considered 
here. At high redshifts, statistical differences in struc- 
ture formation within WMAP3 cosmology when com- 
pared to WMAPl are primarily caused by less small-scale 
power prescribed by the lower erg value (0.9 — > 0.76) and 
scalar spectral index n (1 — 0.96) of primordial den- 
sity perturbations. This manifests in (1) a time delay 

^ See |http: / /Ica.ucsd.edu/software/enzo/ 1 



of ^4 0% of the halo form ation times for a given virial 
mass (jAlvarez et al.l 120061 ) . (2) a corr esponding lower 
halo a bundance for star- forming halos f Gao et al] 120071 : 
IWang e t al. 2008), and (3) stronger clustering of halos 
(jWang et a l. 2008). The initial conditions of this simula- 
tion are well-established by the primordial temperature 
fluctuations in the cosmic microwave background (CMB) 
and big bang nucleosy nthesis (BBN) (|Burles et al.ll200i : 
iHu fc Dodelsonll2002L and references therein). 

We perform two realizations in which we vary the box 
size and random phase to study different scenarios and 
epochs of halo collapse. In the first simulation, we setup 
a cosmological box with 1 comoving Mpc on a side (simu- 
lation A), periodic boundary conditions, and a 128^ top 
grid. The other simulation is similar but with a box 
side of 1.5 comoving Mpc and a different random phase 
(simulation B). We provide a summary of the simula- 
tion parameters in Table [TJ These volumes are adequate 
to study halos of interest because the comoving number 
density of >10* K halos at z = 10 is ~6 Mpc~^ accord- 
ing to an ellipsoidal vari ant of Press-Schechter formalism 
(|Sheth fc Tormenll200l . We use the COSMICS pack- 
age to calculate the initia l conditions'^ at z = 129 (119) 
(|Bertschingeij Il995l 1200 1| ). It calculates the linearized 
evolution of matter fluctuations. We first run a dark 
matter simulation to z = 10 and locate the DM halos 
using the HOP algorithm (|Eisenstein fc Hutlll998l ). We 
identify the first dark matter halo in the simulation that 
has Tvir > 10'' K and generate three levels of refined, 
nested initial conditions with a refinement factor of two 
that are centered around the Lagrangian volume of the 
halo of interest. The nested grids that contain finer grids 
have 8 cells between its boundary and its child grid. Dur- 
ing the simulation, the initial grids retain its position and 
are always refined to its initial resolution or higher. Their 
boundary conditions with each other are treated as any 
other adaptive grid. The finest grid has an equivalent 
resolution of a 1024'^ unigrid and a side length of 250 
(300) comoving kpc. This resolution results in a DM 
particle mass of 30 (101) AIq and an initial gas reso- 
lution of 6.2 (21) A/0. These simulations continue from 
the endpoints of simulations A6 and B6 of Paper I. Table 
[2] lists the parameters of the most massive halo in each 
realization. We evolve the system until the central ob- 
ject has collapsed and reached our resolution limit. If we 
were to follow the simulation to later times and focus on 
subsequently collapsing halos, the nature of the gaseous 
collapses in these halos should be similar because we do 
not consider any non-local feedback processes that affect 
neighboring halos. At redshift 15, the mean separation 
of halos with Tvir > 10"* K is 540 and 910 comoving kpc 
in WMAPl and WMAP3 cosmology, resp ectively, using 
Sheth-Tormen formalism (|Sheth fc Tormen 2002) . Thus 
we argue that the results presented here should be appli- 
cable to all high-redshift protogalactic collapses. 

There are 1.23 x 10* (498^) and 7.40 x 10^ (420^) 
unique cells in the final simulation output of simulations 
A and B, respectively. The finest grid then has a refine- 
ment level of 41 and a spatial resolution of roughly 0.01 
of a solar radius in both simulations. 

Enzo employs a non-equilibrium chemistry model 

* To simplify the discussion, simulation A will always be quoted 
first with the value from simulation B in parentheses. 
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Fig. 1. — An overview of the final state of the collapsing protogalactic gas cloud. Slices of log gas density in cm~^ are shown through the 
densest point in the halo. The first and three rows show simulation A, and the second and fourth rows show simulation B. The columns in 
the top two rows from left to right are slices with a field of view of 10 kpc, 1 kpc, 100 pc, and 1 pc. For the bottom two rows, the fields of 
view are O.Olpc, 20AU, 0.2AU, and 4 Rq. Note that each color scale is logarithmic, spans 5 orders of magnitude, and is unique for every 
length scale. 



()Abel et al.lll997l: lAnninos et al.lll997D . and we consider 
six species in a primordial gas (H, H"*", He, He"*", He"''"'", 
e~). Compton cooling and heating of free electrons from 
the CMB and radiative losses from atomic cooling are 
computed in the optically thin limit. At high densities 
in the halo cores, the baryonic component dominates the 
material. However, the discrete sampling of the DM po- 
tential by particles can become inadequate, and artificial 
heating (cooling) of the baryons (DM) can occur. To 
combat this effect, we smooth the DM particles in cells 



with a width <0.24 (<0.36) comoving pc, which corre- 
sponds to a refinement level of 15. 



3. RESULTS 

In this section, we first describe how the halo collapses 
when it starts to cool through Lya line emission. Then 
we discuss the role of turbulence in the collapse. Last 
we describe the rotational properties and stability of the 
halo and central object. 
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Electron Fraction 



Temperature 





Fig. 2. — Slices of electron fraction (left) and temperature (right) of simulation A (top) and B (bottom). The field of view is 1.5 kpc 
(left panels) and 200 pc (right panels). The color scale is logarithmic for electron fraction and linear for temperature in units of 10^ K. 
Supersonic turbulent shocks are ubiquitous throughout the halos. 



3.1. Halo Collapse 

Beginning at z — 21.1 in simulation A, the progenitor 
of the final halo (Mvi,. = 4.96 x lO^M©) starts to ex- 
perience two major mergers, which continues until z = 
17.2 when Mvir = 2.36 x IO'^Mq. We define Mvir as the 
mass M200 in a sphere that encloses an average DM over- 
density of 200. In simulation B, no recent major merger 
occurs before the cooling gas starts to collapse, but it 
accumulates mass by accretion and minor mergers. 

Mergers disrupt the relaxed state of the progenitor and 
create turbulence as these systems collide and combine. 
Additional turbulence arises during virialization, as dis- 
cussed in Paper I. More small scale density fluctuations 
are thus present in simulation A. These fluctuations pen- 
etrate farther into the potential well in simulation A to 
scales^ of 1 pc, compared to simulation B that contains 
nearly no fluctuations between 1 and 50 pc. This is ap- 
parent in the Z = 1 pc panels of Figure [T] that show the 
density slices at eight length scales covering 11 orders of 
magnitude. At the 10 kpc scale, the filamentary large- 
scale structure is shown, and the protogalactic halo exists 
at the intersection of these filaments. In the next scale, 
we show the protogalactic gas cloud. At the 100 pc scale, 
a thick disc is seen in simulation B. It is nearly edge-on 
and oriented northwest to southeast in this view. In sim- 
ulation B at 1 pc, a bar forms from a rotational secular 
instability that transports angular momentum outwards. 
Similar instabilities exist at radii of 0.2 pc, 2700 AU, 17 
AU, 0.5 R0 in simulation B. Simulation A also under- 
goes a secular bar instability at smaller scales at radii of 
150 AU, 1.3 AU, 0.8 R© but shows a more disorganized 
medium at larger scales. 

^ Note that all masses concerning the collapse are gas mass, not 
total mass. The central regions of r < 10 pc are baryon dominated 



so that M, 

unless otherwise noted. 



one, gas 



' Afenc, tot. AU length scales are in proper units 



The virial temperatures are now > 10* K, and there- 
fore they can efhciently cool by atomic hydrogen tran- 
sitions. The gas fulfills the critical condition for con- 
traction, idyn > icooi, and proceeds to continuously col- 
lapse on approximately a dynamical time. We note that 
this collapse and level of fragmentation are strongly in- 
fluenced by the magnitude of radiative cooling that the 
gas can acheive. Here we present the case in which the 
gas cools without any external radiation backgrounds or 
radiation trapping, which may alter the nature of the 
collapse. 

Figure [2] depicts slices of electron fraction and gas tem- 
perature at scales of 200 and 1500 pc. At the larger scale, 
the gas is heated both in virial shocks at r ^ 600 pc and 
internal turbulent shocks. Gas within the virial radius 
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Fig. 3. — The ratio of the enclosed gas mass and Bonnor-Ebert 
mass (eq. [TJ for the final output (black with circles) and selected 
previous times that are listed in the legend. Simulation A (left) and 
B (right). For values above the horizontal line at Mbuc/Mbe = 1, 
the system is gravitationally unstable. 
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log Radius [pc] '09 ^^nc ["^suJ log Radius [pc] logMgn.IM^J 

Fig. 4. — Mass-weighted radial profiles at various times of (a) gas mass enclosed, (b) number density, (c) mass-weighted temperature, 
and (d) mass-weighted radial velocity for simulation A {left panels) and simulation B {right panels). The quantities in the left and right 
panels are plotted with respect to radius and gas mass enclosed, respectively. In (b), the dashed line with crosses is the dark matter density 
in units of cm"''. In (d), the dashed line with crosses is the negative of the sound speed in the final output. The times in the legends 
correspond to time before the end of the simulation. 
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Fig. 5. — Same as Figure|4]for the inner parsec of simulation A {left panels) and simulation B {right panels). The maximum AMR level 
is listed next to the times in the legend. In simulation B, the local minima in radial velocities at 2 X 10"', 40, 0.3, and 0.01 Mq occur as 
angular momentum is transported outwards in secular bar-like instabilities. 



varies between ^2000 K in cold inflows from filaments 
and up to 30,000 K in turbulent shocks. Electron frac- 
tions increase to up to 0.5% because of coUisional ion- 
izations behind the shocks. The majority of the ioniza- 
tions occur in the turbulent shocks inside rvir where the 
densities are greater and temperatures at the shocks are 
similar to values in the virial shock. In addition, 84% 
of the cooling radiation originates in converging flows 
(V • z; < 0). In the inner 200 pc, turbulent shocks are 
widespread as seen in the temperature variations. How- 
ever these are less pronounced than the one at larger 
radius. In the central 50 pc, the gas becomes nearly 
isothermal despite the low free electron fraction. 

The halo collapses in two stages. We denote the be- 
ginning of the first stage when tdyn > icooi for the first 
time. The second stage begins when the central object 
becomes gravitationally unstable. 

1. Cooling stage — As mass infalls toward the center, 
the increased cooling rate, which is cx nrie until Lya radi- 
ation becomes trapped within the inner condensat ion at 
a density of ~ 5 x 10^ cm'^ (jOh fc Haimanll2002t l. cat- 
alyzes the collapse as atomic line transitions convert ki- 
netic energy to radiation. Here n and rie are the number 
density of baryons and electrons, respectively. Although 
we do not treat the radiative effects of Lya, radiation 



trapping from recombination lines cannot prevent the 
collapse (iRees Sz Ostrikcr 1973) • This first stage starts 
520 (40) kyr before the last output. The inner 100 pc 
have a steady decrease in electron fraction that indicates 
atomic hydrogen cooling is now efficient in this region, 
which can be seen in the 200 pc slices of Figure [2l How- 
ever, only the gas within 1.5 (1.0) pc has idyn ^ ^cooi = 
383 (100) kyr at this epoch. 

2. Gravitationally unstable stage — This starts when 
the ce ntral region b eco mes unstable to gravitational col- 
lapse. lEbertj (|1955l ) and iBonnod ()1955| ) investigated the 
stability of an isothermal sphere with an external pres- 
sure Pext and discovered that the critical mass (BE mass 
hereafter) for gravitational collapse is 

MBE = 1.18^p-y'MQ. (1) 
If we set Pext to the local pressure, then 

Mbe ~ 20r3/^n-i/V"^7^Afo, (2) 

where 7 = 5/3 is the adiabatic index. For both simula- 
tions, this stage occurs between 10 and 100 kyr before 
we end the simulation. We plot the ratio of the enclosed 
gas mass and BE mass in Figure [3] for several epochs in 
the collapse. When the clump becomes gravitationally 
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unstable, the central 3.3 x lO'' (5.5 x 10"^) Mq in the 
central tbe — 5.8 (0.9) pc exceeds the BE mass, and its 
tdyn = 520 (80) kyr. Thus our numerical results agree 
with these analytic expectations. 

We follow the evolution of the accretion and contrac- 
tion until the simulation^ reaches a refinement level of 41 
(41) that corresponds to a resolution of 0.01 (0.014) Rq. 
At this point, the central 4.7 x 10^ (1.0 x 10^) Mq 
are gravitationally unstable and not rotationally sup- 
ported. The central mass is nearly devoid of free elec- 
trons where the electron fraction, ne/n < 10^®, and the 
temperature is ~ 8000 K. It has a radius of 7.9 (1.5) pc. 
The central number density is 5.8 (7.6) x 10^^ cm~^. 
We repeat that this isothermal collapse occurs through 
atomic hydrogen cooling only, but in reality, H2 cooling 
is important even in the presence of a ultraviolet back- 
ground (e.g. iMachacek et al.l l200H f Wise fc Abe]| |2007bt 
10 'Shea fc Normanll2008f) . Thus our results should only 
be considered as a scenario for excellent numerical ex- 
periments of turbulent collapses (see §4.31 for more dis- 
cussion) . 

Next we show the radial profiles of the final and pre- 
ceding outputs in Figures [4] and [5l where we plot (a) en- 
closed gas mass, (b) number density, (c) mass-weighted 
temperature, and (d) mass-weighted radial velocity. Fig- 
ure [4] focuses on length scales greater than 20 AU to 
r > rvir. The halo collapses in a self-similar manner 
with p{r) cx r^'^^/5_ We also overplot the DM density 
in units of mu cm~'^ in the b panels. The DM density 
in simulation A does not flatten as much as simulation 
B with pdm oc r"'*/'^ and r~^/'^, respectively, yet higher 
DM resolution simulations will be needed to address the 
significance of this difference in central slopes. In the 
c panels, ones sees that the entire system is isothermal 
within 10% of 8000 K. In the d panels, the sound speed Cg 
in the final epoch is plotted, and there is a shock where 
Vr > at a mass scale when Mcnc first exceeded Mbe- 
Here Vr is the radial velocity, and Cg is the local sound 
speed. 

Figure [5] shows the data within 1 pc at times 100 
years before the end of the simulation. The self-similar, 
isothermal collapse continues to stellar scales. However, 
the structure in the radial velocity in simulation B ex- 
hibits a strikingly behavior with four repeated minima at 
mass scales 2 x 10**, 10^, 6, and IO'^Mq. We attribute 
this to rotational bar-like instabilities that we discuss 
later in the paper f ji3.6p . 

If we consider Vr constant from the last output, we can 
determine the infall times, which are shown in Figure [6l 
The infall time, tin = r/vr, of the shocked BE mass is 350 
(50) kyr. The infall times approximately follow a broken 
power law, t,n cx Mf„,. Within M^nc ~ O.IM©, /3 « 1/2. 
In the range 0.1 < Mo„c/Mq < 3 x 10^, /3 « 1; above this 
mass interval, the slope of the mass infall times increase 
to P « 3/2. The increased radial velocities when the 
central object becomes gravitationally unstable causes 
the steepening of the slope at ~3 x 10^ M©. 



3.2. Global Disc 

^ We stop the simulation because of ensuing round-off errors 
from a lack of precision. We use 80-bit precision arithmetic for 
positions and time throughout the calculation. 
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Fig. 6. — Radial profiles of gas infall times at the final output. 
To approximate a collapse timescale, the quantities rjv^ (solid) 
and M/ {dM/dt) (dashed) are calculated and plotted here. 

In simulation B, a thick disc with a radius of 50 pc 
and disc scale height of ~10 pc forms that is pressure 
supported and only partially rotationally supported. The 
circular velocities within this disc achieve only a third of 
Keplerian velocities. The lack of full rotational support 
and large scale height suggests that a central collapse 
occurs before any fragmentation in this large-scale disc 
is possible. In contrast, we see a disorganized, turbulent 
medium and no large scale disc formation in simulation 
A. 

3.3. Turbulence 

iKolmogorovl (|1941l ) described a theory of the basic be- 
havior of incompressible turbulence that is driven on a 
large scale and forms eddies at that scale. These eddies 
then interact to form smaller eddies and transfer some of 
their energy to smaller scales. This cascade continues un- 
til energy is dissipated through viscosity. In supersonic 
turbulence, most of the turbulent energy is dissipated 
through shock waves, which minimizes the local nature 
of cascades found in incompressible turbulence. 

In Paper I, we found that turbulence is stirred during 
virialization. When radiative cooling is efficient, the gas 
cannot virialize by gaining thermal energy and must in- 
crease its kinetic energy in order to reach equilibrium, 
which it achieves by radial infall and turbulent motions. 

In addition to virial turbulence generation, mergers stir 
turbulence. Here the largest driving scale will be approx- 
imately the scale of the merging objects, and the turbu- 
lent cascade starts from that length scale. Additional 
driving may come from Kelvin-Helmholtz instabilities of 
a multi-ph ase gas as the mergers occur (, Takiz awa„ 2 5.) . 
iTakizawal considered mergers of galaxy clusters, however 
his work may still apply to the formation of protogalac- 
tic halos since similar temperature contrasts exist in this 
regime of mergers. As the lighter halo falls into the mas- 
sive halo, a bow shock and small-scale eddies from the 
Kelvin-Helmholtz instability form between the two inter- 
acting objects. At later times, a dense, cool core remains 
in the substructure of the lesser halo. The instabilities 
grow and destroy the baryonic substructure, and the gas 
mixes with the existing gas of the massive halo and be- 
comes turbulent. 
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Fig. 7. — The turbulent Mach number, 

^rms I Cs , for the final 

output (black with diamonds ) and selected previous times that are 
listed in the legend. Simulation A {left) and B {right). 

To quantify aspects of this turbulence, we inspect the 
turbulent Mach number. 



M = 



dP jkT 



(3) 



Cs " dp nmn' 

Here P is pressure, Vrms is the 3D velocity dispersion, 
and 7 is the adiabatic index that we set to 5/3. We 
evaluate Vrms with respect to the mean velocity of each 
spherical shell. Radial profiles of Ai are shown in Fig- 
ure [T] Before the core becomes gravitationally unstable, 
the turbulence is subsonic within the virial shock. Af- 
ter the core becomes gravitationally unstable, the turbu- 
lent Mach number rises to 2-4. The collapse produces 
turbulence on a timescale that is faster than it can be 
dissipated. 

The turbulence that exists before the initial collapse 
may impact the nature of the central object. In simula- 
tion A, the core initially has Al « 1, and this results in a 
central object with 4.7 x IO^Mq and a radius of 7.9 pc. 
The core in simulation B has Ai « 0.2, and the central 
object is about five times less massive and smaller, which 
corresponds to a free-fall time approximately five times 
shorter as well. 

3.4. Spin Parameter Evolution 

During the hierarchical buildup of structure, tidal 
forces from neighboring structures impart angular mo- 
mentum to a halo, particularly w hen its radiu s is max- 
imal at the turn-around time (jHovld 1 19491 : iPeeblesI 
119691 ). However in recent years, several groups have 
recognized that the mergers may impart a consid- 
erable ^fracti^n^of^jmgular momentu m to th e sys- 
tem (iSteinmet z fc Bartclmann 199^ iGardneH I20i31l : 
iVitvitska et al.i,2002; Mailer et al. 2002 ). Over many re- 
alizations of mergers, the net angular momentum change 
would be zero. In reality, an angular momentum residual 
remains after the last major merger occurs because there 
are too few events to cancel the randomization of halo 
spin. Although each halo has unique rotational proper- 
ties, it is useful to define a dimensionless spin parameter 



where G is the gravitational constant and L, E, and M 
are the angular momentum, energy, and mass of the 
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Fig. 8.— Spin parameter, A = \L\^/\E\/GM^/^ , evolution of 
the main halo in the simulation, (left) simulation A. (right) sim- 
ulation B. The dashed and solid lines are the interpolated values 
for the DM and baryonic spin parameter. The squares and cir- 
cles correspond to the actual measurements from the DM and gas 
data, respectively. The horizontal dashed line at A = 0.04 marks 
the mean cosmological spin parameter. In simulation A, two ma- 
jor mergers causes the large increase beginning at z 21 in the 
hashed region. The oscillations occur as the merging halos orbit 
each other until they virialize. 

objec t, that measu res the rigid body rotation of the 
halo (iPeebled I1971D . In Figure [51 we display the time 
evolution of A of the DM and baryons in our simu- 
lations and mark the occurrence of the ma jor merger 
in simulation A. lEisenstein fc Loebl (|1995b( l (preceded 
bv lBarnes fc Efstathioulll987D calculated that the mean 
spin parameter, (A) « 0.04, is weakly dependent on ob- 
ject mass and cosmological model, and this value is also 
marked in Figure [H Also A weakly depends on its merger 
history, where (A) increases during mergers and slowly 
dissipates afterwards. Most of the angular momentum 
is acquired from steady minor mergers and accretion be- 
cause major mergers only happen rarely (usually only 
once per logarithmic mass interval). In 96% of mergers, 
the majority of the internal spin originates from the or- 
bital energy of the infalling halo (jHetznecker fc Burkertl 
[200l . 

At z ~ 22 in simulation A, the spin parameter A = 0.06 
before the last major merger. Then the spin parameter 
increases by a factor of 3 during its major merger because 
of the system being far from dynamical equilibrium. The 
system becomes virialized after approximately a dynam- 
ical time, and the spin parameter stabilizes at A ~ 0.03 
and proceeds to decrease with time until A = 0.022 at 
the time of collapse. The above evolution of A agrees 
with the findings of Hetznecker & Burkert. Simulation 
B describes a halo that does not undergo a recent major 
merger, and its final A — 0.013. 

Both halos have less angular momentum than (A) when 
the cooling gas collapses. The probability distribution of 
A can be described with the log-normal function 



p{X)dX = 



1 



exp 



ln'(A/Ao) 



2crx 



where Ap = 042 ± 0.006 and ax = 0.5 ± 0.04 (e.g. 
iBullock et al.l |2001[ ) . From the cumulative probability 
function resulting from equation ([5]), 89% (99%) of the 
cosmological sample of halos have l arger spin parame- 
ters than the halos described here. lEisenstein fc Loea 
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Fig. 9. — Secular instability e-folding timcscalc in units ofaf/u 
as a function of t = T/\W\ and a = (t//2)^/2 (gq^ [gj ^t t 
< 0.1375, the system is stable to all perturbations. Above t = 
0.27, the system is dynamically unstable, and this timcscalc is not 
applicable. 



(|1995af) demonstrated that halos with low spin parame- 
ters are candidates for BH formation and quasar seeds. 
However they argue that the angular momentum needs 
to be at least an order of magnitude lower than the mean. 
Next we present further evidence that reveals a gaseous 
collapse is possible with not too atypical spin parameters. 

3.5. Instability of Maclaurin Spheroids 

The dynamics of rotating systems is a classic topic in 
astrophysics (see EFE §§1-6). These self-gravitating sys- 
tems are susceptible to two types of instabilities. Secular 
instability occurs when small dissipative forces, e.g. vis- 
cosity, amplify perturbations to become unstable in an 
otherwise stable inviscid configuration. Dynamical (also 
referred to as ordinary) instability results when some os- 
cillatory mode exponentially grows with time, regardless 
of any dissipative forces. Here we concentrate on Maclau- 
rin spheroids relevant for a uniform body rotating with a 
fixed angular velocity. Maclaurin spheroids are a special 
case of Jacobi ellipsoids that are axisymmetric. The on- 
set of the m = 2 bar-like instability in gaseous Maclaurin 
spheroids happens for a given eccentricity, 



1 



1/2 



> 



0.8127 
0.9527 



(secular) 
(dynamical) 



(6) 



where 03 and ai are the principle axes with 03 < oi (EFE 
§33). Eccentricity is related to the ratio, t — T/\W\, of 
rotational kinetic energy to gravitational potential by 



t 



2)-3(t-=-l)"^(sm-'e)->l, (7) 



and the secular and dynamical inst abilities happen at 
t = ( 0.1375,0.27), respectively (e.g. lOstriker fc PeeblesI 
[l973h . 

When t is larger than 0.1375 but smaller than 0.27, 
both the Maclaurin spheroid and Jacobi ellipsoid are 
perfectly stable against small perturbations in the in- 
viscid case. For a given e, the Jacobi configuration has 
a lower total energy than its Maclaurin counterpart and 



is therefore a preferred state. Here any dissipative force 
induces a secular bar-like instability. The system slowly 
and monotonically deforms through a series of Riemann 
S-type ellipsoids until its final state of a Jacobi ellipsoid 
with an equal angular momentum (iPress fc T eukolskvl 
[l973) and lower angular velocity (EFE §32) as specific 
angular momentum is transported outward. The insta- 
bility grows on an e-folding timescale 



(t>al/v, 



(8) 



where is a constant of proportionality that asymptotes 
Sii t ~ 0.1375, decays to zero &l t = 0.27, and is plotted 
i n Figure [9l (EFE §37). Here v is the kinematic viscosity. 

IChristodoulou et all ()1995al lbl) generalized the formu- 
lations for bar-like instabilities to account for self-gravity. 
In addition, they consider different geometries, differ- 
ential rotation, and non-uniform density distributions. 
They devised a new stability criterion 



where f2 is the rotation frequency, 
"(l-e2)i/2 




■ sm 



(9) 



(10) 



is the Jeans frequency in the radial direction for a 
Maclaurin spheroid, and 



/ 



1 



1 



(11) 



accounts for differing geometries^ with / = 2/3 for a 
sphere and / = 1 for a disc. Secular and dynamical insta- 
bilities for Maclaurin spheroids occur above a = (0.228, 
0.341), respectively, for / = 1. 

From N-body simulations of disc galaxies, 
lOstriker fc Peeblei ()1973D found that a massive dark 
halo with comparable mass to the disc could suppress 
secular instabilities. In the case of a gaseous collapse to 
a SMBH however, the baryonic component dominates 
over the dark matter component in the central 10 pc. 
Secular instabilities cannot be prevented through this 
process, which we demonstrate next. 

3.6. Rotational Instabilities 

In the Z = 1 pc panel of simulation B in Figure [l] it is 
apparent a bar-like instability exists in the gravitation- 
ally unstable central object. Figure [TOl shows the insta- 
bility criterion a (eq. ^ against enclosed gas mass. Here 
we transform the velocities to align the z-axis with the 
baryonic angular momentum vector of the entire halo. 
We use the tangential velocities to calculate the rota- 
tional kinetic energy T. The shape parameter / = 2/3 
(0.89) for simulation A (B). 

As discussed before, Maclaurin spheroids are subject 
to secular m = 2 bar-like instabilities when a > 0.228. 
In simulation A, the central object becomes unstable at 
three approximate mass scales, 6.7 x 10""*, 1.0, and 110 
Mq that correspond to radii of 0.75 R©, 1.3 AU, and 150 
AU, respectively. The enclosed mass ratios of the recur- 
ring instabilities, i.e. M^/Mi+i, are 1500:1, 110:1, and 



See 

tries 



IChristodoulou et al.l l|1995bl ') for more generalized geome- 
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Fig. 10. — Rotational instability parameter a = fT/2\W\ for 
the thick disc with r ~ 50 pc in simulations A (black solid line) and 
B {red dashed line). The shaded areas show the standard deviation 
when varying the center on the 100 densest points. For a > 0.22 
denoted by the horizontal line, a secular instability occurs in the 
disc and leads to bar formation. In simulation A, instabilities occur 
at mass scales of 100, 0.1, and 10~*Mq. In simulation B, the same 
happens at 2 X 10^, 2 X 10"', 10^, 6, and IO'^Mq. We also mark 
a = 0.341 where a rotating system becomes dynamically unstable. 
Only simulation A at 0.1 Mq experiences a dynamical instability. 

1400:1, starting at the smallest mass scale. The instabil- 
ity at 0.075 Mq (r = 0.13 AU) is dynamically unstable 
with a peaking at 0.55. In simulation B, instabilities oc- 
cur at 5.3 X 10"-^, 7.0, 1.2 x 10^, and 2.0 x W^Mq at 
radii of 0.49 R©, 17 AU, 2700 AU, and 0.18 pc. The 
enclosed mass ratios of these instabilities, are 13,000:1, 
170:1, 17:1, and 85:1. 

It is interesting to note that the innermost instability 
in both simulations becomes dynamical (a > 0.341), and 
a continues to increase rapidly toward the center. How- 
ever these features should be taken with caution since 
it occurs near our resolution limit, where the particular 
location used as the center will influence the rotational 
energy one would calculate. To evaulate the sensitivity 
in choosing a center, we performed the same analysis but 
varying the center over the 100 most densest cells in the 
simulation. We plot the standard deviation of a as the 
shaded area in Figure [TOl Inside an enclosed mass of 
3 x 10"'' M©, it is ^0.05 but diminishes to less than 0.01 
outside O.IMq. 

The e-folding time of secular instabilities r is propor- 
tional to of (see eq. Hence small-scale instabilities 
collapse on a faster timescale than its parent, large-scale 
bar instability. Turbulent viscosity is the main dissipa- 
tive force that drives the instability, r is inversely pro- 
portional to the viscosity. This further shortens r be- 
cause supersonic turbulence is maintained to the smallest 
scales. 

3.7. Rotational Properties 

During the collapse of the gas in our simulations, ro- 
tational support never impedes the collapse. In Figures 
[TT] and [121 we show (a) coherent rotational velocity di- 
vided by Keplerian velocity Vkep — \/ GM jr, (b) rota- 
tional velocity, (c) specific angular momentum, and (d) 
rotational velocity divided by the sound speed. We com- 
pute the rotational velocities around the center of mass 



of a sphere with radius of 100 cell widths of the finest 
AMR level, centered on the densest point. We note that 
the rotational velocity L/r plotted here is different than 
organized rotation, i.e. a disc. The radial proflles only 
sample gas in spherical shells, whose angular momentum 
vectors are not necessarily parallel. 

1. Simulation A— At r > 1 AU (Menc = l-M©), the 
typical rotational speed is two or three times lower than 
the Keplerian velocity, which is required for rotational 
support. At r = 0.1 AU (Mc„c = 0.07 Mq), the infaU 
becomes marginally rotationally supported, i.e. L/r ~ 
Vkep- The radial velocities react by slowing from 15 km 
s~' to below 5 km s~'. However this rotational support 
does not continue to the center. Rotational speeds are 
only r-^Q.bvkep within 0.1 AU (Monc = IMq). 

2. Simulation B — This collapse exhibits four min- 
ima in radial velocity that are caused by rotational bar- 
like instabilities. After such an instability occurs, the 
radial velocities increase because of angular momentum 
being transported outwards. As the rotational veloci- 
ties decrease, this instigates another secular instability, 
which repeats causing a cascade of the instability. The 
increased infall velocity and associated decrease in rota- 
tional velocities (i.e. the dips in Figures [SJi and [l2tl) 
depict this behavior. At the final output, the infalling 
material exhibits no rotational support at all radii simi- 
lar to simulation A at r > 1 AU. 

We interpret the inner points where L/r/ukep fluctua- 
tions greatly increases above unity with caution because 
of the nature of choosing a center in a turbulent medium, 
i.e. when turbulent velocities dominate over rotational 
ones. If the central sphere is smaller than a radius where 
the turbulent velocities average to zero, we introduce er- 
rors into the angular momentum profiles by sampling 
the turbulent gas incompletely. In the 6-panels of Fig- 
ure llli one sees that specific angular momentum inside 
Monc < IO^Mq decreases over time and is transported 
outwards in the collapse. 

With a not too atypical spin parameter, the thick disc 
with r ^ 50 pc is not rotationally supported. In simula- 
tion A, a global disc does not exist at all. We attribute 
this behavior to the nature of angular momentum trans- 
port in a turbulent medium. Even with a higher spin pa- 
rameter, we do not expect a disc to fragment before the 
central collapse of gas with low specific angular momen- 
tum and short dynamical times. This low specific angu- 
lar momentum material collapses to small radii without 
fragmentation so that a central dense object forms with 
a mass of ~ IO^Mq or 2% of the halo gas mass. After the 
initial collapse, the thick disc may become rotationally 
supported as more high angular momentum gas infalls. 

4. DISCUSSION 

In our cosmological simulations, we find that a 
~10^Mq dense object forms in the center of a metal-free 
protogalactic halo that cools by atomic hydrogen cooling. 
Although we have neglected some important processes, 
such as H2 chemistry, star and BH formation and feed- 
back, our results show that angular momentum transport 
at both small and large scales in the form of preferential 
segregation and rotational instabilities, respectively, lead 
to the formation of a dense, massive object with r < 5 
pc. This initial central collapse should precede any frag- 
mentation of a global disc. 
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4.1. Angular Momentum Transport 

Collapsing turbulent clouds, whether cosmological or 
galactic in nature, are ubiquitous in the universe. In this 
paper, we focus on the details of the turbulent collapse 
of a proto-galactic halo. Angular momentum transport 
plays a key role in such events, e.g., determining the char- 
acteristics of the central object(s). However there exists 
the "angular momentum problem" , where many orders 
of magnitu de of angular momentum must be shed (see §6 
in lLarsonir2003l from the initial molecular cloud to form 
a central star, star cluster, or BH. In our simulations, 
there is a clear s cenario in which the inside-out collapse 
(|Shu et al.l I1987D proceeds even if the initial turbulent 
cloud was rotating. We see three major elements affect- 
ing angular momentum transport during the collapse. 

1. Angular momentum distribution — In cosmological 
halos, there is a universal distribution of angular momen- 
tum 

M(<j) = M™.-^, M>1, (12) 

Ja+J 

that measures the mass with a specific angular momen- 
tum less than j (iBuUock et al.l |2001[ ) . This function 
is fitted with two parameters, /i and joi where fi con- 
trols the flattening of the power law at high angular mo- 
menta, and Jo determines at which j this transition oc- 



curs. iBullock et al.l also find that more mass resides in 
the tails of the distribution, especially at small j, when 
compared to a system in solid body rotation. Thus all 
halos have some intrinsic amount of gas with small j . If 
this distribu tion is maintained during the collapse (e.g. 
lMestellfl963| ). such gas can collapse to some small radius, 
'''min > j/vkep, without bccoming rotationally supported, 
which leads to the next element of discussion - angular 
momentum segregation. 

2. Segregation in a turbulent medium — In Paper I, 
we determined that most of the gas becomes superson- 
ically turbulent as a result of virialization. Therefore 
let us theorize how angular momentum transport hap- 
pens during the transition from being pressure supported 
to rapidly cooling and collapsing. First consider a tur- 
bulent uniform-density gas cloud, where parcels of gas 
at a specific radius can have many different values of j. 
This differs from the organized rotation of a disc. If we 
start with such an initial configuration, how does angu- 
lar momentum transport occur during the collapse? Gas 
with small (high) j will preferentially migrate to small 
(large) radii, following turbulent flow li nes. In an axi- 
symmetric system, th e Rayleigh criterion (jRavleighlll92Cll : 
lChandrasekhMl[l96l requires that the specific angular 
momentum must be a monotonically increasing function 
with respect to radius. The gas with the lowest j pro- 
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gressively piles up in the center of DM potential wells 
until tcooi < tdyn when it can catastrophically cool and 
collapse. Such low j gas may originate in lower mass 
progenitors because the gas resided in shallow potential 
wells (i.e. low mass halos) that led to smaller turbu- 
lent and thermal velocities. We argue that this effect is 
intimately linked to the gas acting to achieve virial equi- 
librium at all stages during the collapse (see Paper I). 
Furthermore, the system becomes unstable to turbulence 
as the material segregates. This onset of turbulence can 
be delayed if viscosity is large enough so that Reynolds 
numbers are below the order of 10^ or 10^. However 
there are many modes of instability if the Rayleigh crite- 
rion is not met, and even gas with a low Reynolds num- 
ber will eventually become fully turbulent on a timescale 
that is chaotic, dep ending on t he initial perturbatio n and 
Reynolds number (piulll992t iMoehhs et al.ll2004D . We 
note that a more comprehensiv e approach would con sider 
the Solberg-H0iland criterion (jEndal fc Sofialll978D that 
generalizes this to include partial rotational and pressure 
support in a disc. 

3. Bar-like rotational instabilities — After sufficient 
amounts of gas have migrated to small radii because of 
angular momentum segregation, this gas increases its ro- 
tational velocity as it conserves angular momentum. Gas 
with similar angular momentum now obtains some orga- 
nized rotational velocity. As the rotational energy in- 
creases, some shells may become rotationally unstable 
(T/|W| > 0.14) in a secular m — 2 mode. In the case of 
a collapsing gas cloud, turbulent viscosity provides the 
dissipative force that drives the secular instability. The 
system then deforms into a bar-like object, where the gas 
with large j moves to larger radius and gas with small j 
can infall to even smaller radii. 

The combination of these three processes alleviates the 
"angular momentum problem" of inside-out collapses. 
Such a scenario of angular momentum transport dur- 
ing a self-similar collapse may be widely applicable in 
cosmological collapse problems. 

4.2. Secular Instability Cascade 

Our simulations follow the self-similar collapse of pro- 
togalactic halos over 14 orders of magnitude in length. 
We find that a cascade of three (four) bar-like instabili- 
ties occur during the latter stages of the collapse. The ra- 
tios of mass enclosed in each successive instability varies 
from 10 to 10,000 in our simulations. As a consequence 
of these instabilities, the collapse of the densest point 
never halts because of rotational support. Instead the 
gas becomes rotationally unstable when it gains suffi- 
cient rotational energy. The lowest j gas then falls to 
smaller radius and may become unstable yet again. This 
sequence could repeat itself several times. In addition, 
we find that rotational instabilities are possible without 
a global disc as in simulation A. 

This is the "bars within bars" scenario originally pro- 
posed to fuel active galactic nuclei through dynamical ro- 
tational bar-like in stabilities (jShlosman et al.l[l989lll99dl : 
iHeller et al.l[200l . It was then adapted for tunneling 
enough gas into pre-ga lactic (Af ~ IO^Mq) SMBHs by 
iBegelman et al.l ( 20061 ). in whose framework the angu- 
lar momentum of the disc, where the instability occurs, 
depends on the spin parameter of the halo (see also 



iMo et aLlll998D . Thus the amount of gas available for ac- 
cretion onto the central SMBH also depends on the spin 
parameter. Dynamical instabilities require 45% more ro- 
tational energilto_occur than secular ones. In the frame- 
work of IBegelman et aP . only requiring secular instabil- 
ities may result in a larger fraction of halos forming a 
pre-galactic SMBH because of the log-normal distribu- 
tion of spin parameters (eq. [5]). Nevertheless, we do 
not advocate our simulations as evidence of pre-galactic 
SMBH formation because we have neglected many im- 
portant processes related to H2 cooling and primordial 
star formation that we detail briefiy in the next section. 

4.3. Applicability 
4.3.1. Limitations of Current Approach, 

Our results depict the importance of turbulence, accre- 
tion, and the hydrogen cooling in the initial collapse of 
these halos. However we are missing some essential pro- 
cesses, such as H2 chemistry, primordial and Population 
II stellar formation and feedback, SMBH formation and 
feedback, and metal transport and cooling. It was our 
intention to study only the hydrogen and helium cooling 
case first and gradually introduce other processes at a 
later time to investigate the magnitude and characteris- 
tics of their effects, which we will present in later papers. 

Gas becomes optically thick to Lya radiation above 
column densities of ~10^'^ cm~^, and Lya radiation 
trapping becomes important a bove a density of ~5 x 
10 cm~ 3 ()0h fc Haimanll2002f) . We continue to use op- 
tically thin cooling rates above this density. Thus we 
overestimate the cooling within 0.03 pc. As a conse- 
quence, we do not suggest that these simulated objects 
ever form in nature. However this scenario poses an ex- 
cellent numerical experiment of turbulent collapse, which 
should be common in galaxy formation, where turbu- 
lence is generated during virialization, and star forma- 
tion within turbulent molecular clouds. 

4.3.2. Desired Improvements 

Clearly local dwarf spheroidals contain stars with 
ages consistent with formation at very high red- 
shifts (iFerrara fc Tolstovll2000l : lTolstov et al.|[200a l2003t 
iHelmi et al.l 120061 ) . To develop a model that desires to 
fit galaxy luminosity functions down to the faintest ob- 
served galaxies one may need a star formation and feed- 
back model that follows molecular clouds as small as one 
thousand solar masses in order to allow for the dominant 
mode of star formation observed locally. It should be al- 
ready technologically feasible with current cosmological 
hydrodynamical models to simulate these galaxies one 
star at a time. 

Correct initial conditions for early galaxy formation 
require prior star and BH formation and feedback. 
The typically adopted conditions for phenomenologi- 
cal star formation are velocity convergence, a criti- 
cal overdensity, ^dvn > ^ cooii and being Jeans unsta- 
ble (jCen fc Ostrikerl Il992| ). Phenomenological primor- 
dial star formation is possibl e if we include tw o addi- 
tional conditions as utilized in I Abel e t al.' (2007*). First, 
the H2 fraction must exceed 10~^ (Abel et al. 2002), 
and second, the mctallicity of the gas must not ex- 
ceed some " critical metalhcity" of 10~'^ - 10~^ of the 
solar value (|Bromm et al.l 120011: ISchneider et all 120061 : 
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ISmith k Sigurdssoiill2007l:lJappseii et al.ll2007allH). From 



prior studies (e.g.lAbel et al .1120021: iBromm fc Loebl [20?)3l: 
lO'Shea et al] 120051 : iGreif &: Bromml I2006D . we expect 
these stars to form to in halos that can support II2 cool- 
ing and ones embedded in relic H II regions. The 
Lyman-Werner radiation from m assive stars can dis- 
sociate H2 from la rge distances (jPekel fc ReesI llQSTt 
iHaiman et al.ll2000D. suppress sta r formation in lowe r 
mass halos (jMachacek etall [2OOII: iWise fc Abel [2001 . 
and should be considered to accurately model future star 
formation. 

BH formation in the death of some primordial stars 
can also have a profound effect on surrounding struc- 
ture formation as it accretes infalling matter during later 
mergers. In principle, one should include feedback of 
seed BHs from primordial stars with masses outside of 
the range between 140 and 260 solar masses. Also it 
is possible to phenomenologically model SMBH forma- 
tion in a similar manner as the stellar case. If the pro- 
togalactic collapse occurs faster than stellar formation 
timescale of a massive star, a SMBH may form inside this 
region. Using the stellar formation conditions plus this 
condition and allowing the particle to accrete (i .e. sink 
particles; iBate et al.iri995l: iKrumholz et al.ll2004) , proto- 
galactic collapses can be followed in cosmological hydro- 
dynamic simulations (Clark ct al. 2007). These sink par- 
ticles should regulate the accretion with an appropriate 
subgrid model. Important processes include an appro- 
priate accr etion rate (e.g. Eddin gton or Bondi-Hoyle), 
turbulence (|Krumholz et al.ll2006[ ). rotational support of 
the infalling gas, and a viscosity timescale for accretion 
discs. 

For small galaxies , radiative transfer effects can have 
a great impact ( e.g. iHaehneltlllQQSl: iWhalen et al.l l2004l: 
iKitavama et al.l 12004': 'Alva rez et all l2006 f) and should 
not be neglected. Ionization front instabilities in these 
galaxies create cometary small-scale structure and shad- 
owing effects as a result from the explicit treatment 
of three-dimensional radiation hydrodynamics. Stellar 
feedback can have both a positive and negative impact 
on subsequent star formation. Some examples of positive 
feedba ck include enhanced H2 formation i n relic H II re- 
gions (|Ferraral [l998l : lO'Shea et all 120051 : iJohnson et ail 



20071) and dust and i netal-line cooli r ig (iGloved 120031 : 
Schneider et~al] 120061: iJappsen et all l2007al) . Nega- 



tive feedbac k may occur f rom b aryonic expulsion from 
host halos (IWhalen et~ ar 2004!; IKitavama et al.l 12004 
lYoshida et al.l 120071: I Abel et al.i 120071) and halo pho- 
toeva poration ( Susa fc Umemural 120061 : IWhalen et all 
[2OO8I) . 

The promising approach of lGnedin fc Abel (|2001l ) has 
recently been implemented and coupled with t he AMR 
hydrodynamic code ART (|Gnedin et al.l l2008l). Also 
the t echnique of adaptive rav tracing (jAbel fc Wandelg 
I2OO2D has been implemented into Enzo and used to 
study the out flows and ionizing radiation from a pri- 
mordial star (l Abel et all 120071 ). This method has 
also been independently impleme nted into Enzo by 
iRazoumov fc Sommer-LarsenI (I2006D. Finally as used in 
many stellar form ation routines ( Cen fc Ostrikerl 119921 : 
iTassis et al.ll2003l ), we hope to include thermal and ra- 
diative feedback from Population II stars in future stud- 



ies. 



5. CONCLUSIONS 

We have simulated the hydrodynamics and collapse of 
a protogalactic gas cloud in two cosmology AMR realiza- 
tions. Our focus on the hydrodynamics presents a basis 
for future studies that consider stellar and BH feedback. 
In the idealized case presented, we find a central dense 
object forms on the order of lO^M© and r ;< 5 pc. This 
central object is not rotationally supported and does not 
fragment in our simulations. However our results do not 
dismiss disc formation in protogalaxies because rotation- 
ally supported disc formation may begin after the initial 
central collapse. Disc formation may be sensitively af- 
fected by feedback from the central object. 

These simulations highlight the relevance of secular 
bar-like instabilities in galaxy formation and turbulent 
collapses. Similar bar structures are witnessed in pri- 
mordial star formation simulations. As low angular mo- 
mentum infalls, it gains rotational energy as it conserves 
angular momentum. This induces an m = 2, bar- like in- 
stability that transports angular momentum outwards, 
and the self-similar collapse can proceed without becom- 
ing rotationally supported and exhibits a density profile 
p (X r~^^/^. This process repeats itself as material infalls 
to small scales that is indicative of the "bars within bars" 
scenario. We see three and four occurrences of embedded 
secular instabilities in the two realizations studied here. 

We also find that supersonic turbulence influences the 
collapse by providing a channel for the gas to preferen- 
tially segregate according to its specific angular momen- 
tum. The low angular momentum material sinks to the 
center and provides the material necessary for a central 
collapse. Here the possibilities of a c entral object in- 
clude a direct collapse into a SMB H (e.g. Bromm fc Loebl 
l2003h . a starbur st (e.g.jC lark et al.l i2007f ). or a combina- 
tion of both (e.g. lSilk fc R ces 199^ All of these cases are 
viable in the early universe, and the occurrence of these 
cases depends on the merger history, local abundances 
in the halo, and the existence of a seed BH. Moreover, 
star formation should occur whether a central BH exists 
or not. Perhaps the frequency of these different proto- 
galactic outcomes may be traced with either 3D numer- 
ical simulations that consider star and SMBH formation 
and feedback along with metal transport or Monte Carlo 
merger trees that trace Pop HI star formation, metallic- 
ities, and BHs. We will attempt the former approach in 
future studies to investigate protogalactic formation in 
more realistic detail. 
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